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Abstract 

We investigate the effects of low-lying fermion eigenmodes on the QCD partition 
function in the e-regime. The fermion determinant is approximated by a truncated 
product of low-lying eigenvalues of the overlap-Dirac operator. With two flavors of 
dynamical quarks, we observe that the lattice results for the lowest eigenvalue distri- 
bution, eigenvalue sum rules and partition function reproduce the analytic predictions 
made by Leutwyler and Smilga, which strongly depend on the topological charge of the 
background gauge configuration. The value of chiral condensate extracted from these 
measurements are consistent with each other. For one dynamical quark flavor, on the 
other hand, we find an apparent disagreement among different determinations of the 
chiral condensate, which may suggest the failure of the e-expansion in the absence of 
massless Nambu-Goldstone boson. 
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§1. Introduction 



Chiral perturbation theory (ChPT) is an effective field theory to describe the dynamics 
of Quantum Chromodynamics (QCD) at low energies (<C A ~ 1 GeV), where the Nambu- 
Goldstone pion excitations dominate the dynamics while other particles are too heavy to be 
excited. In the infinite volume, ChPT provides a method to express low energy amplitudes 
of pion as an expansion in terms of pion mass squared m 2 and its momentum squared p 2 . 
Thus it enables us to calculate low energy pion amplitudes in a systematic manner. When 
the system is put in a finite volume, for which the pion Compton wavelength (~ 2rr/m n ) 
is much larger than the linear extent L of the space-time, i.e. m n L <C 1, while L is kept 
large enough compared to the QCD scale 1//1q C d, the low energy effective theory can still 
be constructed as an expansion in terms of small e 2 ~ m w /A ~ p 2 /A 2 , which is known as 
the e-expansion. 1 ) In this setup, the so-called e-regime of ChPT, the chiral symmetry is 
not spontaneously broken and one must explicitly integrate over different vacua in the path 
integral, which leads to characteristic behavior of the partition function and other physical 
quantities. In particular, they strongly depend on the gauge field topology (or the existence 
of fermion zero modes). 2 -* 

Lattice QCD simulation is well suited to the study of finite volume physics. Since the 
chiral symmetry plays an essential role in the e-regime, one should use lattice fermion for- 
mulations that respect chiral symmetry. Such fermion formulation has become available 
relatively recently, i.e. the overlap fermion 3 )' 4 ) and the domain- wall fermion. 5 )' 6 ) These 
fermion formulations satisfy the Ginsparg- Wilson relation, 7 ) with which the fermion action 
can be shown to have exact chiral symmetry at finite lattice spacings. 8 ) Quenched lattice 
simulations have so far been done in the e-regime to calculate the chiral condensate 9 ^ and 
other low energy constants. 10 H 3 ) It has also been shown that the distribution of low- lying 
eigenvalues of the lattice overlap Dirac operator agrees well with expectations from the chiral 
Random Matrix Theory (RMT) 14 )~ 16 ) (for a review of RMT, see 17 - 1 ), which is equivalent to 
the chiral Lagrangian at the leading order of the e-expansion. 

In this work we extend these previous lattice studies to the case of non-zero dynamical 
fermion flavors. In particular we calculate the topological susceptibility and partition func- 
tion with one and two flavors of dynamical quarks, and investigate their dependence on the 
quark mass. The predictions from ChPT in the e-regime are available from the work by 
Leutwyler and Smilga. 2 ) Due to the presence of fermion zero-modes the partition function is 
drastically different for different topological sectors of background gauge field. We therefore 
calculate them for each topological sector identified by the number of fermion zero-modes. 
We also investigate the sum rules of Dirac operator eigenvalues (the so-called Leutwyler- 
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Smilga sum rules 2 -*) derived from the quark mass dependence of the partition function. For 



these quantities we found that the lattice data with one or two fermion flavors are well 
described by the known analytic formulas. 

We also extend the comparison of low-lying eigenvalue distributions with the RMT results 
to the case of non-zero number of flavors. In RMT there is an universality among the 
number of flavor Nf and topological charge v. Namely the prediction depends only on the 
combination Nf + which should be observed in the lattice data. 

In order to investigate the effect of dynamical fermions, we employ an approximation for 
the fermion determinant. Namely, we approximate the fermion determinant by a product of 
low-lying eigenvalues of the overlap Dirac operator. The effects of higher fermion modes are 
neglected, which we call the truncated determinant approximation. Since the eigenvalues 
much higher than the quark mass and the QCD scale Aqcd should be irrelevant to low 
energy physics, this approximation is expected to be effective for the study of the e-regime. 
The higher eigenmodes are sensitive to lattice artifacts and their main effect for low energy 
physics appears through the renormalization of parameters in the theory, i.e. in this case 
the gauge coupling constant and quark masses. As we neglect such effects, we effectively 
choose a different renormalization scheme. To what extent this approximation works can be 
tested by varying the cutoff for the Dirac operator eigenvalues. 

This paper is organized as follows. First, in Section |21 we briefly review the analytic 
expectations from ChPT following the seminal paper by Leutwyler and Smilga. 2 - ) We then 
explain our calculation methods and describe the truncated determinant approximation in 
some detail in Sections El and |U respectively. Numerical results are presented in Section El 
An earlier presentation of this work is found in. 18 ) 



using the chiral Lagrangian in the e-regime. 

At the leading order in the e-expansion the kinetic term in the chiral Lagrangian is 
suppressed, because the momentum excitation energy is large in a small volume even for the 
unit momentum 2n /L. The partition function is then given by 



§2. Leutwyler-Smilga's analytic predictions 



In this section we briefly review the analytic predictions made by Leutwyler and Smilga 2 -* 



Z = exp [EVReime" 



)] 



(2-1) 



for N- 



f - 



1 and 




(2-2) 
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for Nf > 2. Here, E represents the chiral condensate and V is the space-time volume. 
They appear in the combination mEV together with the quark mass m. The integral in 
()2-2|) runs over compact SU(iV/) elements Uq corresponding to the zero momentum mode 
wave function, and d[i(Uo) is its Haar measure. This integral is necessary because the chiral 
symmetry is not spontaneously broken in the finite volume. For Nf = 1 (Eq. ()2TJl ). on the 
other hand, this degree of freedom does not remain due to the U(l) chiral anomaly. The 
parameter 9 is the CP angle of the QCD vacuum. 

The partition function in each topological sector labeled by the topological charge v is 
obtained by Fourier transforming ()2-l|) and ()2-2|) . They are 

Z v = I v (x) for N f = 1, (2-3) 

Z v = I 2 (x) - I v+x {x)I v - X {x) for N f = 2. (2-4) 

I v {x) is the modified Bessel function and its argument x is x = mEV . Near the massless 
limit it behaves as 

The probability to find a configuration with a given topological charge v is given by 
Z v jZ, where Z = ^2 V Z V is e x and I 1 (2x)/x for Nf = 1 and 2, respectively. The topological 
susceptibility (z/ 2 ) /V is then obtained through the definition (z/ 2 ) = ^ v y 1 Z v jZ as 



mE for N f = 1, (2-5) 

for Nf = 2. (2-6) 



V 

(is 2 ) mEI 2 {mEV) 



V 2h{mEV) 

Asymptotically, it behaves as (u 2 ) /V oc m Nf for i C 1 and (u 2 ) /V = mE/Nf for x ^> 1. 

The above results for the partition function should agree with the partition function of 
the underlying theory, QCD, in the appropriate limit. The QCD partition function for a 
given topological charge v is written as 



m 



N f \v\ 



[[dG}e- s ^f[(\\ n \ 2 + m 2 ) N f. (2-7) 



The path integral over the gauge field configuration is done for a fixed topology v with a 
weight given by the gauge action Sq- The fermion determinant is written in the form of a 
product of the eigenvalues A n of the Dirac operator. Non-zero eigenvalues appear with their 
complex conjugate and the index n in (J2-7j) runs over eigenvalues with positive imaginary 
part only. The prime on the product indicates that it excludes the zero modes, which are 
factored out as m N f^. Defining the expectation value in the massless limit as ((■ • • )}„, the 
above expression is rewritten as 

m- N fMZJm) //A A m 2 N N 



lirn^o m- N fMZJm) \\ 1A V |A 



n ^ ■ c- 
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In the effective theory the same quantity is given by 



Mi(§)\(*) 



for Nt = 1, 



(2-9) 



[v + l)\\u- 1 ! 



2|i/| 



Z„(a;) 



for N f = 2. (2-10) 



By taking a derivative with respect to m 2 for both QCD and the effective theory and equating 
them at m = 0, we arrive at the so-called Leutwyler-Smilga sum rules. From the first 
derivative we obtain 

lit \\ 

i \\ i 



Y< (\\ n \zvY/j v mf+w\y (2 ' n) 

which we call the sum rule I in this paper. The prime on the sum implies that it does not 
include the zero modes. For non-degenerate flavors with masses and rrij we can extract 
two sum rules from second derivatives corresponding to mjm 2 and mf 



E 



E 



\\n\ZVf 



E 



1 



lQ((N f + \u\y 



(\\ n \£Vy // 16(N f +\v 
which we call the sum rules Ha and lib, respectively. Their difference gives 

1 \\ 1 



l)(JV> + M+2) 



(2-12) 



(2-13) 



E 



(\x n \Evy 



16(JV / +|z/|)((JV / + |i/|) 2 -l) : 



(2-14) 



(sum rule lie). Note that there is a universality among the number of flavors and topological 
charge, i.e. they appear only in the combination Nf + \v\. 

These sum rules fj211jl - (|2-14)l do not make sense as they stand, because the left hand 
side is ultraviolet divergent while the right hand side is a definite number. In fact, since the 
eigenvalue density p(A) = h (J2 n $(^n ~ A)) increases as ~ V^A 3 for large A, the sum rule I is 
quadratically divergent. The sum rules Ha and lib have quartic divergences, which cancel in 
the rule lie and only a logarithmic divergence remains. In spite of these divergences, the sum 
rules fl2-llj) - ()2-14|) are consistent, because they only affect the corrections of order 1/V or 
higher, and thus can be removed by first taking a limit V — > oo. In the following numerical 
work we do not perform this extrapolation, which is very expensive, but consider differences 
between different topological sectors. Since the eigenvalue density becomes independent of 
topology for large A, the ultraviolet divergence cancels in such differences. 
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Topological charge \u\ 


1 2 


3 


4 


5 


6 


#config 


169 290 149 


68 


20 


4 


3 



Table I. The number of quenched gauge configurations for each topological sector. 



§3. Lattice details 

The chiral symmetry is essential for the lattice study in the e-regime. We use the Neu- 
berger (or the so-called overlap) Dirac operator 3 )' 4 ) which respects the Ginsparg- Wilson 
relation, and thus the chiral symmetry is preserved at finite lattice spacing. 

The overlap Dirac operator D is defined as 

D = -{l + l5 sgn(aH w )], (3-1) 
a 



with 

Hw = 75 D w - 

a 



(3-2) 



D w = \Y, + V„) - aV^V^} , (3-3) 

and a — a/(l + s). Here Dw is the conventional Wilson-Dirac operator and V (V*) denotes 
a forward (backward) gauge-covariant differential operator on the lattice. The parameter s 
is introduced to optimize the negative mass term given in the kernel Hw- In this work we 
choose s = 0.6 at (3 = 5.85. The sign function in ()3T|) is approximated using the Chebyshev 
polynomial after subtracting out several lowest eigenmodes of Hw- The order of polynomial 
is optimized for each gauge configuration to ensure the accuracy of 10~ n for the sign function, 
and it is typically around 100-200 when we subtract 14 lowest-lying eigenvalues of Hw- 

We work on quenched gauge configurations generated with the plaquette gauge action at 
/3 = 5.85 on a 10 4 lattice. The lattice spacing determined through the Sommer scale ro is 
0.123 fm. The topological charge for these gauge configurations is determined by the number 
of zero modes of the overlap-Dirac operator ()3Tj) . The number of gauge configuration for 
each topological sector is given in Table |U 

For the overlap-Dirac operator ()3Tjl we calculate 50 lowest eigenvalues and their eigen- 
vectors for each gauge configuration. We utilize the numerical package ARPACK, 19 - 1 which 
implements the implicitly restarted Arnordi method to compute the eigenvalues. Instead of 
treating the operator D as it is, we consider a chirally projected operator D ± = P±DP±, 
where P± = (1 ±7s)/2. The size of the matrix is then a factor of 2 smaller and the numerical 
calculation becomes about 2 times faster. The eigenvalues of D can be obtained from the 
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eigenvalues of P±DP± } using the property that the chirally projected operator gives a real 
part of the eigenvalues of the original operator. Since the eigenvalues of D lie on a circle 
satisfying |1 — A| 2 = 1, we can construct a pair of eigenvalues of D, namely A and A*. The 
corresponding eigenvectors u of D can also be reconstructed from the eigenvectors = P±u 
of the chirally projected operator using a formula 

u = u + = u . 3-4 

1mA 1mA 

In order to identify the number of left-handed and right-handed zero modes, we have to 
compute the eigenvalues of both P + DP + and P_DP_. We carry out such calculation until 
we achieve enough precision to identify pairing non-zero eigenvalues of both operators. 20 - ) The 
number of zero modes is then determined unambiguously, and we continue the calculation 
of the rest of the eigenvalues on the operator for which we do not find the zero modes. 

§4. Truncated determinant approximation 

In our calculation the effects of dynamical fermions are incorporated in the partition 
function by including a product of eigenvalues of Dirac operator in the pure gauge path 
integral. 



m 



N f \v 



J [dUU~ s -H(Xl + mTf =m*W (j[(Xl + rnYf\ (44) 



where A n = |A n |y 1 — (am) 2 /4. The eigenvalues of the overlap-Dirac operator appear in 
complex conjugate pairs (A, A*) except for the zero-modes as in the continuum case. The 
notation (•••)^ denotes an expectation value evaluated on quenched gauge configurations 
with a fixed topological charge v. Namely, we generate the gauge configurations with a weight 
determined by the pure gauge action and reweight them with the product of eigenvalues. 

This method requires a calculation of all eigenvalues of the overlap-Dirac operator for 
each configuration, but it is not feasible unless the lattice size is extremely small (~ 4 4 ). 
Instead, we approximate the fermion determinant by a truncated product of eigenvalues in 
(j4-lj) . The upper limit of the eigenvalue plays a role of an ultraviolet cutoff for the fermionic 
degrees of freedom. Qualitatively, this procedure is justified, since the physical effects of the 
fermion determinant are expected to come from the low-lying eigenmodes, while the higher 
eigenmodes should be irrelevant for low energy physics. As in the usual renormalization 
program, cutoff dependence of physical quantity can be absorbed in the coupling constants 
and anomalous dimensions up to lattice artifacts. In our case the bulk of the effects of higher 
fermion modes is neglected, and thus the relation between the lattice bare coupling and the 
lattice spacing is mostly the same as in the quenched approximation. 
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The role of the cutoff in the fermion determinant can be rephrased by defining a new 
Dirac operator Dim (LM stands for low-lying modes) as 



D LM = 75/ A cu tofr) (4-2) 
with H = 75D and a function f(x, A cu toff) satisfying a condition 

{X (\x\ <C Acutoff): 

U 1 cutoff, (43) 
Acutoff (pi ^> Acutoff)- 

The truncation corresponds to a non-smooth function, which transfers from x to a constant 
Acutoff at Acutoff- One can also interpolate the two regions smoothly using an analytic func- 
tion, e.g. /(x, Acutoff) = Acutoff fanh(x/A cu toff)- For such analytic functions the locality and 
unitarity of the new Dirac operator Dlm can be proved. 2l > Furthermore, the operator -Dlm 
satisfies a modified Ginsparg- Wilson relation 

7sD L m + -Dlm75 = aD^ 5 D LM = aD LM j 5 D, (4-4) 

as far as f(x, A cu t ff) is written in terms of a polynomial. It implies that the fermion action 
constructed with the truncated operator -Dlm is invariant under the chiral transformation 
5ip = 75(1 — ^P-)il> and 5tp = ip(l — ^)75, just like the original overlap-Dirac operator. The 
truncated operator -Dlm m ay, therefore, be used as an alternative definition of the Dirac 
operator with exact chiral symmetry. It should be noted, however, that our choice for the 
function f(x, A cu toff) (hard cutoff) has to be understood as an approximation as it is not a 
smooth analytic function. 

On a 10 4 lattice at (3 = 5.85, we calculated 50 smallest eigenvalues of the chirally projected 
overlap-Dirac operator. Among them, the maximum eigenvalue corresponds to the physical 
scale ~ 1200 MeV, which is large enough as a cutoff for the study of the low energy physics 
of interest. We define the truncated determinant as 

A^ut 

det(D> + m) = Yl C^l + ™ 2 ) (4-5) 

n=l 

with a fixed number A^ cu t of eigenvalues. Figure Q shows a correlation between the truncated 
determinant with N cut = 10 and that with N cnt = 50 for topological sectors v — 0, 1, and 2. 
The quark mass is set to zero, and the zero-modes are factored out for non-trivial topological 
sectors. We observe that they are strongly correlated up to an order of magnitude varia- 
tions, while the truncated determinant itself varies as much as 8 orders of magnitude when 
measured on quenched gauge configurations. It suggests that the truncated determinant 
provides a good approximation of the full determinant already at -/V cu t = 10. 
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1e-21 1e-20 1 e-19 1e-18 1e-17 1e-16 1 e-15 1e-14 1e-13 

Truncated determinant with 10 eigenvalues 
Fig. 1. Correlation between the truncated determinant (|4-5|) with iV cut = 10 (horizontal axis) 
and that with N cnt = 50 (vertical). Each point corresponds to a gauge configuration with 
topological charge (pluses), 1 (crosses), and 2 (stars). The quark mass is set equal to zero. 
Zero-modes for nontrivial topological sectors are not included in the product. 
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Fig. 2. Relative weight due to the truncated determinant for the massless overlap-Dirac operator 
as a function of N cnt for five representative gauge configurations with v = 0. 



One may also see how much the relative size of the truncated determinant depends on 
the value of cutoff N cut . In Figure we plot the truncated determinant for the massless 
overlap-Dirac operator as a function of N cnt for five representative gauge configurations, for 
which we calculated 300 smallest eigenvalues. The value is normalized by an average over 
these five configurations at a given value of N cnt , and thus the plot shows a relative size 
of the weight factor due to the approximate determinant. We find that the relative weight 
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varies rapidly for small iV cu t (~ 10), while it becomes roughly constant above iV cu t > 50. 
It supports our expectation that the truncated determinant is a good approximation of the 
full determinant for N cnt ~ 50. In the following studies we also test how much our results 
depend on the parameter N cut for each quantity we measure. 

The number of necessary eigenvalues to be included in the truncated determinant is 
expected to be unchanged even at smaller lattice spacings as long as the physical volume 
is kept fixed. This is because the the eigenvalue density p(A) depends on the physical 
parameters £ and V but not on the lattice cutoff 1/a (up to lattice artifacts). It depends on 
the physical volume linearly, and therefore the number of eigenvalues to be calculated would 
become prohibitively large, e.g. about 14 times larger on a (2 fm) 3 x(4 fm) lattice, which is 
a typical lattice size used for hadron spectrum calculations. 

The truncated determinant was previously introduced by Duncan, Eichten and Thacker 
to develop an efficient algorithm for light dynamical quarks with the Wilson fermion. 22 ) They 
also proposed a method to include the higher fermion modes to make the algorithm exact 
using the multiboson technique of Luscher. 23 ) A similar method can also be applied in our 
case, but the effective bosonic action becomes non-local with the overlap-Dirac operator and 
the simple heat-bath updation would be impractical. 

An immediate problem of the reweighting with the approximate determinant is that the 
overlap of the unquenched vacuum with the quenched vacuum could be tiny and thus the 
Monte Carlo sampling becomes inefficient especially for small quark masses we are interested 
in. If this happens, only a few configurations give a dominant contribution to the Monte Carlo 
average and others are suppressed by the reweighting factor representing the approximate 
determinant. We may quantify this effect by considering an effective number of statistics 
iVeff as 

N«= (4-6) 
maxj(det-D ) 

(0 



Here, i and j label gauge configurations and iV is the number of the configurations, det D 



(j) 

denotes the truncated determinant for the z-th gauge configuration, and maXj(det-D ) is 
its maximum value in the gauge ensemble. Figure El shows the effective number of statistics 
(normalized by the quenched statistics) as a function of quark mass. The effective number 
of statistics decreases for smaller quark masses as expected. For a given topological sector, it 
goes down to ~ 0.05 (0.02) for Nf = 1 (2), when the quark mass is decreased to am ~ 0.01, 
which corresponds to m ~ 26 MeV. If we combine all topological charges, we loose another 
factor of ~ 5, because only the configurations with zero topological charge contribute for 
small quark masses. Therefore, the number of configurations we are using (Table |l} are not 
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ALL TOP.SECTOR 

v=0 
|v|=1 
|v|=2 



0.1 



0.01 



0.001 
0.0001 



0.001 



0.01 
am/(1+s) 



0.1 



ALL TOP.SECTOR 

v=0 

|v|=1 
|v|=2 



N, = 2 



0.1 



0.01 



0.001 1 — — — — 

0.0001 0.001 0.01 0.1 1 

am/(1+s) 

Fig. 3. Effective number of statistics N c g normalized by the quenched statistics N for Nf = 1 
(top panel) and Nf = 2 (bottom panel). It is calculated for each topological sector (squares, 
crosses and bursts for \u\ = 0, 1, and 2, respectively) and for all topological sectors (pluses). 



large enough for precision study of unquenched lattices. 



§5. Numerical results 



In this section we present our numerical results for the eigenvalues of the overlap-Dirac 
operator. An eigenvalue A of the overlap-Dirac operator lies on a circle defined by |1— a\\ = 1. 
In the continuum limit a — > they become pure imaginary, corresponding to the eigenvalues 
of the continuum Dirac operator. At the finite lattice spacings, on the other hand, there is 
a freedom to define the eigenvalue projected onto the imaginary axis. We simply take an 
absolute value of the eigenvalue of the overlap-Dirac operator. The ambiguity in the choice 
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5 10 15 20 25 30 35 

AXV 

Fig. 4. Eigenvalue distribution p{\)/ S of the overlap-Dirac operator as a function of XEV. Dis- 
tributions are shown for topological sectors \u\ = (solid), 1 (dashed) and 2 (dotted). The 
asymptotic behavior ~ A 3 in the bulk region is shown by a solid curve, while the fits with RMT 
are drawn by dashed curves. 

of the projection leads to a systematic uncertainty of order (a\) 2 . In the following discussion 
we misuse A for the meaning of |A|. iV cu t is 50 unless otherwise stated. 

5.1. Eigenvalue distribution 

First of all, we show the distribution of the low-lying eigenvalues in Figure 0] The plot 
shows the eigenvalue density p(X) = J2 n 8(^n — A)) normalized by the chiral condensate 
E. This normalization is determined for each topological sector by comparing an average of 
the lowest eigenvalue with its RMT expectation, as discussed below. 

As expected from the dimensional analysis, the spectral density behaves as ~ A 3 in the 
large XEV region. This behavior is independent of the gauge field topology, and a fit with 
Co + c^XEV) 3 for the entire gauge field ensemble is shown by a solid curve (co = 0.55, 
C3 = 1.15 x 10~ 4 ). Below XEV ~ 10, on the other hand, one can see a difference among 
different topological sectors. The curves show a fit with the expected functional form from 
the chiral RMT. 

In Figure we compare the distribution of the lowest non-zero eigenvalue with the pre- 
diction from RMT. 24 ) Unlike the previous works, 14 -* we also obtain the distributions for 
non-zero number of flavors, Nf = 1 and 2. For the quenched case, our results confirm the 
previous works that found an agreement with the RMT predictions. 

The unquenched distributions are obtained by the reweighting with the truncated deter- 
minant. The lowest eigenvalue distribution is probably one of the most difficult quantities 
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Fig. 5. Lowest non-zero eigenvalue distribution on gauge configurations with a topological charge 
v = (left), 1 (middle) and 2 (right). The results are shown for Nf = (top), 1 (middle), 
and 2 (bottom). The horizontal axis is normalized such that the average agrees with the RMT 
expectation. 



to be measured using the reweighting method. For instance, for v = almost all eigenvalues 
lie in the region XEV < 5 on the quenched lattice, while only about a half of eigenvalues 
are expected to fall in that region for Nf = 2 as shown in the plot by a dotted curve. The 
other half in the region XEV > 5 is not sampled well with the reweighting. 

In Table ITT1 we summarize the average of the lowest non-zero eigenvalue Ai predicted by 
RMT and our results. For the RMT an expectation value of \\SV is listed, while the lattice 
data are bare values aX\. From these results we can extract the value of chiral condensate 
E for each \u\ and Nf. For the lattice spacing we use the quenched value a = 0.123 fm even 
for unquenched data, assuming that the heavy quark potential is not much affected by the 
truncated determinant. The chiral condensate £ obtained in this work corresponds to the 
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RMT (AiZV) 


Lattice (aAi) 


Z" 1 / 3 [MeV] 




n 

u 


1 77 


n n9RQf'i i\ 




Nt— o 


1 


3 1 1 


n n486fi 2^ 


2^4(2^ 




2 


4.34 


0.0701(20) 


251(2) 




n 

U 


^11 


u.uouy ^oy j 


z^i j 


l\f— 1 


i 
± 






9^8 1"}"! 




2 


5.53 


0.100(6) 


242(5) 







4.34 


0.067(12) 


254(15) 


iV>=2 


1 


5.53 


0.0939(42) 


247(04) 




2 


6.69 


0.112(14) 


248(10) 



Table II. Expectation value of the lowest non-zero eigenvalue predicted by RMT X\UV and the 
lattice results aAi. The chiral condensate E is obtained with an input for lattice spacing a = 
0.123 fm at (3 = 5.85. 



bare lattice operator ipifj and is not renormalized. To convert it to the continuum definition, 
such as the MS scheme, we need a matching factor. 

The agreement of S given in Table ITT1 among different topological sectors is encouraging. 
The values with non-zero flavors are also in reasonable agreement, though they are not 
necessarily agree since the underlying theory is different. 

5.2. Partition function 

Employing the truncated determinant approximation, we calculate the unquenched par- 
tition function at a fixed topology Z v (12 -7|) as a function of the quark mass. 

The number of eigenvalues included in the truncated determinant has to be equal for 
different topological charges in order to make the mass dimension consistent for differ- 
ent topological sectors. This is relevant for the partition function, because we consider 
Z v jZ, and Z in the denominator is a sum over all topological charges. Therefore, we 
slightly modify the truncated determinant as rw v v\ is even, 

or as ml"' ]Jn=i~^ +1 ^ 2 fin + rn 2 ) x ■\J^ 2 Ncut _Q 1/ ^ 1 y 2 + 777,2 wnen M is °dd. 

In Figure we plot Z v jZ as a function of mEV for Nj = 1 (top panel) and 2 (bottom 
panel). We find a good agreement with the expectation from the chiral Lagrangian in the 
small quark mass region ()2-3|) and ()2-4|) . A fit in the range dm < 0.03 with a free parameter 
E yields A7 1 / 3 = 214(6) MeV and 248(6) MeV for N f = 1 and 2, respectively. These results 
are quite insensitive to the value of N cnt . They vary from 212(4) MeV to 214(6) MeV by 
changing N cut from 10 to 50 for N f = 1. The variation is from 234(7) MeV to 248(6) MeV 
for N f = 2. 
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Fig. 6. Partition function for a given topological charge as a function of the quark mass. The 
results are shown for Nj = 1 (top) and 2 (bottom). The curves represent the analytic prediction 
in the e-regime. 



5.3. Eigenvalue sum rules 

Here we investigate the sum rules (|2Tljl - (|2T4jl of the Dirac operator eigenvalues. 

As discussed in Section E] the sums Y^' n V-^n an d J2 n are divergent in the ultraviolet 
regime. This behavior is explicitly shown in Figure [7| for Nf = 1. The left panels are 
divergent sums as a function of A cut ZV with A cut the largest eigenvalue included in the sum. 
The sum rule I (top left) does not seem to show the expected quadratic divergence but looks 
more like a linear divergence. But, it is in fact consistent with the integral J Xcut dXp(X)/X 2 
with the form p(A) = p^\X) + p^A 3 suggested in Section I5~T1 Here, p^A 3 represents the 
bulk distribution of eigenvalues, which is independent of the gauge field topology, while the 
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Fig. 7. Sum rules I (|2-lll) and lib ()2-13|) for Nf = 1. The ultraviolet divergence is shown on the 
left panels as a function of X cn t£V with A cut the largest eigenvalue included in the sum. The 
right panels show a saturation for the difference between different topological sectors. 



pi (A) term shows the oscillating behavior depending on the topological charge. For large 
XSV, p^\X) approaches to a constant p(°\ which is independent of the topology. In the 
integral J Xcut dXp(X)/X 2 , we observe that the curvature of —p^/X cut cancels that of p^A^. 
in the region plotted in Figure [7| 

Similar plots are shown for Nf = 2 in Figure |H1 The first order (sum rule I) and the second 
order (sum rules Ha, lib and lie) are plotted. The sum rule I is quadratically divergent as 
in the Nf = 1 case; the divergence of the sum rule lie ()2T4j) is mild as it is only logarithmic. 
For the other second order sum rules (Ha and lib), the divergence is quartic, since it has the 
form ( f Acut dXp(X)/X 2 ) 2 . 

In order to subtract these ultraviolet divergences, we consider the differences of the 
eigenvalue sum rules among different topological sectors. This is motivated by an expectation 
that the bulk distribution p^A 3 is independent of the gauge field topology as we can see in 
Figure HJ On the other hand, the term p^ carries the information of the topology and shows 
the oscillating behavior expected from RMT. For the sum rule I (j2Tlj) . to be explicit, the 
divergent piece are common for all topological sectors and cancel in the differences. In fact, 
the lattice data for the differences [v — 0) — {y = 1), {y — 0) — [y — 2) and [y = 1) — [y = 2) 
are almost independent of the cutoff A cut as shown in Figure (top right panel). The same 



16 



(v=0)-(v=1) 



(v=0)-(v=2) 



(v=1)-(v=2) 



sum rule I 



sum rule Ha 



sum rule lib 



sum rule lie 



V=0 ' 

v=1 > 
v=2 



v=0 > 
v=1 > 
0.06 - v=2 : 

0.05 - 

0.04 

0.03 

0.02 

0.01 





0.05 




0.045 


> 






0.04 




0.035 


E 


0.03 








0.025 


> 


0.02 








0.015 


E 


0.01 


en 


0.005 



v=0 > 

V=1 H 

v=2 



v=0 ^ 
v=1 H 
v=2 



: 1 ■ j. i i .!. ■ ■ ' 



15 20 25 30 35 



| 0.04 




1 



15 20 25 30 35 






0.02 








0.045 




04 




0.035 




0.03 


> 






0.025 


O 


0.02 


E 


0.015 


CO 






0.01 




0.005 









0.03 


> 






0.025 


o 


0.02 


E 




=> 
CO 






0.015 


> 






0.01 


o 
E 


0.005 








5 10 15 202530 35 5 10 15 20 25 30 35 



(v=0)-(v=1) 





5 10 15 20 25 30 35 
(v=1)-(v=2) 



5 10 15 202530 35 5 10 15 20 25 30 35 5 10 15 20 25 30 35 



(v=0)-(v=1) 




(v=1)-(v=2) 



15 20 25 30 35 



| | | | i ^ 

5 10 15 202530 35 5 10 15 20 25 30 35 5 10 15 



20 25 30 35 



(v=0)-(v=1) 



(v=0)-(v=2) 



(v=1)-(v=2) 



C 0.01 
o 

E 

0.005 



10 15 20 25 30 35 



5 10 15202530 35 5 10 15 20 25 30 35 5 10 15 20 25 30 35 

K,J» K„£V 



Fig. 8. Sum rules I (t2~TTT) . Ila (t2~12ll . lib (THTil) and lie (t2~TH) for N f = 2. The ultraviolet 
divergence is shown on the left panels as a function of X cn tUV with A cu t the largest eigenvalue 
included in the sum. The right panels show a saturation for the difference between different 
topological sectors. 
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Fig. 9. Expectation value of the topological charge squared (u 2 ) as a function of mSV. Results 
for Nf = 1 (pluses) and 2 (crosses) are shown. The dashed line and the dotted curve are fit 
results with ()2-5|) and ()2-6|) . respectively. 



is true for the sum rule lie as shown in Figure |H1 (4th row, right). 

The cancellation of the ultraviolet divergence is not expected for the sum rules Ha ()2T2)1 
and lib (|2T3j) . This is because they include a square of the integral (/ cut d\p{\)/\ 2 ) 2 and 
it contains a cross term pi, > (A) x p^\ 2 ut , which is divergent but still topology-dependent 
through pi\\). The lattice data support this expectation as shown in Figure (bottom 
right panel) and in Figure |H1 (2nd and 3rd rows, right). 

We use the sum rules I and lie to extract £ from the comparison of lattice data with 
the analytic formulae. For Nf = 1, only the sum rule I can be used, and we obtain E 1 ^ 3 
as 239(9), 238(7) and 233(9) MeV for the differences [y = 0) - (v = 1), (v = 0) - (y = 2) 
and (u — 1) — (v — 2), respectively. Similarly, for Nf = 2 we obtain 268(17), 260(13) and 
241(17) MeV from the sum rule I, and 257(12), 254(11) and 237(10) MeV from the sum rule 
He. 

5.4. Topological susceptibility 

The mass dependence of the topological susceptibility (v 2 )/V in the e-regime is analyti- 
cally known as ()2-5|) and ()2-6|) . Lattice calculation of this quantity is straight-forward from 
the partition function discussed in Section 15.21 In Figure El we plot (u 2 ) as a function of 
m£V. The value of E extracted by fitting in the region am < 0.03 is Z 11 / 3 = 216(9) MeV for 
Nf = 1, and 257(18) MeV for Nf = 2. These results are quite stable under the change of iV cut . 
The variation for iV cut = 10-50 is 215(7)-216(9) MeV for N f = 1 and 246(27)-257(18) MeV 
for N f = 2. 
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The clear mass dependence of the topological susceptibility as observed in Figure El has 
not been obtained before in the dynamical fermion simulations, 25 ^ except an exploratory work 
by Kovacs, 26 - 1 which utilized the same reweighting technique as employed in our work. The 
reason is partly that the quark mass is not small enough to realize the e-regime. A recent work 
with slightly smaller quark masses shows an indication of the suppression of the topological 
susceptibility. 27 - ) Also, the (improved) Wilson fermion formulation employed in the previous 
dynamical simulations does not respect the chiral symmetry, and it is questionable if the 
expected quark mass dependence is obtained at finite lattice spacings. For the improved 
staggered dynamical quarks, with which one can reach the small quark mass region, there 
is an indication that the expected quark mass dependence is reproduced in the continuum 
limit 28 ) 

5.5. Comparison 

In Table HTTl we summarize the values of chiral condensate E extracted from our numerical 
data. We set a = 0.123 fm from the Sommer scale ro at (3 = 5.85. 

First of all, the quenched (Nf = 0) results obtained from the lowest eigenvalue distribution 
are consistent among different topological sectors, and in perfect agreement with a previous 
work by Bietenholz et a/., 15 ) who used the same lattice parameters, i.e. (3 = 5.85 and 10 4 
lattice, and extracted E from the lowest eigenvalue distribution as E 1 / 3 ~ 253 MeV. In 13 ) 
we investigated the meson correlators in the e-regime at the same (3 value but on a larger 
lattice 10 3 x20. The quenched chiral condensate is extracted from the scalar and pseudo- 
scalar meson correlators as E 1 ^ 3 ~ 257(14) MeV, which is also consistent with the present 
work. 

For Nf = 1, we find an agreement between the results from the lowest eigenvalue distri- 
bution and those from the sum rule I. The results are around E 1 ^ 3 ~ 240 MeV, which is in 
the same ballpark with the quenched results. On the other hand, the fits for the partition 
function and the topological susceptibility yield substantially lower values, 214(6) MeV and 
216(9) MeV, respectively. At first sight, this disagreement seems puzzling, because the eigen- 
value sum rules are derived through the derivatives of the partition function with respect 
to m 2 . The main difference is that the partition function and the topological susceptibility 
involve the relative weight among different topological sectors. Namely, the leading mass de- 
pendence of the partition function Z v oc m'"', which is most relevant to the relative weight, 
is factored out before the sum rules are derived. The eigenvalue sum rules represent the 
non-leading mass dependence. Therefore, the disagreement implies that the relative weight 
among different topological sectors is no well described by (|2-3|) . 

In contrast to the Nf = 1 case, the results for E from several different methods nicely 
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Table III. Summary of the results for chiral condensate I7 1 / 3 . Results from the lowest eigenvalue 
distribution, eigenvalue sum rules I and lie, partition function, and topological susceptibility 
are listed for Nf = 0,1 and 2. The second column indicates the topological charge of the gauge 
configuration used for the measurement. 

agree with each other for Nf = 2. Probably, the extraction from the lowest eigenvalue 
distribution contains substantial systematic error due to the sampling problem, as discussed 
in Section 15.11 The sum rules are less problematic, because they involve the inverse of 
the low-lying eigenvalues and its most important contribution comes from lower side of the 
distribution, which is well sampled with the reweighting method. 
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§6. Conclusions 



The effect of the fermion determinant to the QCD vacuum is substantial for small enough 
quark masses. In the language of chiral perturbation theory, the constant modes dominate 
the low energy dynamics, when the system is confined in a finite volume. The quark mass 
dependence of the partition function can be analytically derived using the e-expansion in 
the chiral perturbation theory. The analytic formulae by Leutwyler and Smilga show strong 
dependence on the topological charge and the number of dynamical quark flavors. In this 
numerical work we explicitly tested these analytic expectations using lattice QCD simula- 
tions. 

The fermion determinant is approximated by a truncated product of low-lying eigen- 
values of the overlap-Dirac operator. The eigenvalues are truncated at around 1200 MeV. 
Intuitively, this approximation should be effective as far as the low energy dynamics is con- 
cerned, and in fact the relative weight among gauge configurations is roughly unchanged 
when we further increase the cutoff for the eigenvalue. 

From our lattice calculations of (i) lowest-lying eigenvalue distribution, (ii) eigenvalue 
sum rules, (iii) quark mass dependence of the partition function, and (iv) topological sus- 
ceptibility, we found good agreement with the analytic predictions. For these quantities, 
the only free parameter in the chiral perturbation theory (at the leading order) is the chiral 
condensate E, which we can extract from the lattice data. The extraction of E from the 
methods (i)-(iv) provides a highly non-trivial cross-check of the theory (or the lattice cal- 
culation), and our results at Nf = 2 are all consistent with each other. A problem remains 
at Nf = 1, for which the extraction with a fixed topological charge ((i) and (ii)) and the fits 
to all topological charges ((iii) and (iv)) are mutually inconsistent. It may signal a break- 
down of the e-expansion at Nf — 1, for which no Nambu-Goldstone massless pion exists and 
therefore the finite momentum excitation is not large compared to the ground state energy. 

Although the reweighting method with the truncated determinant provides a good quali- 
tative approximation, its limitation is also apparent. First of all, the systematic uncertainty 
due to the truncation is hardly estimated. For the quantities (i)-(iv) we explicitly checked 
that the results are unchanged by varying the number of eigenvalues from 10 to 50, but it 
does not guarantee that they remain unchanged until the limit of the exact determinant is 
reached. This problem could be solved either by defining an effective Dirac operator which 
includes the low-lying eigenvalues only or by incorporating the effects of the rest of the eigen- 
modes using stochastic methods. Second, the Monte Carlo sampling with the pure gauge 
action becomes ineffective when quark mass becomes small. Near the massless limit, we 
loose nearly 95% (98%) of the statistics for Nf = 1 (2) by the suppression factor due to the 



21 



fermion determinant. This problem can only be cured by including the fermion determinant 
into the Boltzmann weight when the gauge configuration is generated. 
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